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ABSTRACT 

High resolution dark matter only simulations provide a realistic and fully general 
means to study the theoretical predictions of cosmological structure formation models 
for gravitational lensing. Due to the finite number of particles, the density field only 
becomes smooth on scales beyond a few times the local mean interparticle separation. 
This introduces noise on the gravitational lensing properties such as the surface mass 
density, the defiection angles and the magnification. At some small-scale mass limit, 
the noise due to the discreteness of the N-body simulation becomes comparable to 
the effects of physical substructures. We present analytic expressions to quantify the 
Poisson noise and study its scaling with the particle number of the simulation and the 
Lagrangian smoothing size. We use the Phoenix set of simulations, currently the largest 
available dark matter simulations of clusters to study the effect of limited numerical 
resolution and the gravitational strong lensing effects of substructure. We quantify 
the smallest resolved substructure, in the sense that the effect of the substructure on 
any strong lensing property is significant compared to the noise, and we find that the 
result is roughly independent of the strong lensing property. A simple scaling relates 
the smallest resolved substructures in a simulation with the resolution of the N-body 
simulation. 

Key words: Gravitational lensing - methods: numerical - dark matter - galaxies: 
clusters: general 



1 INTRODUCTION 

The cold dark matter (CDM) paradigm makes two major 
predictions on the mass structure of dark matter haloes. 
Dark matter haloes at all scales are expected to have a uni- 
versal mass density profile and to be populated by a large 
number of mass substructures (Gao et al. (2004); Diemand 
ct ai. (zuU4); bprmgei ct ai. (z005); Gao et al. (2012)). Grav- 
itational lensing provides a unique tool to probe the mass 
distribution of galaxies and galaxy clusters and therefore to 
test this model. At galaxy cluster scale, a number of large 
gravitational lensing surveys (e.g. CLASH, Postman et al. 
(2012); SLOAN, Oguri et ai. (201^)) have accurately deter- 
mined the shape of the mean mass density profile. Using a 
combination of weak and strong lensing data for example 
Nowmmi ot (2009) and TTmofsn et al. (2011) have mea- 
sured the mass distribution of galaxy clusters from kpc to 
Mpc scales. In general, both in simulations and observations, 
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the central profile of clusters shows a significant amount of 
scatter between different observations (Sand et al. (2002, 
2004); Newman et al. (2011)) as well as between different 
simulations (Gao et al. (2012)); this scatter might be ex- 
plained in the context of a hierarchical structure formation 
model where clusters form late and therefore are not fully 
relaxed objects. 

Combining weak and strong lensing is especially impor- 
tant to measure the cluster concentration, defined as the ra- 
tio of the virial radius to the radius where the density profile 
has an isothermal slope. To date, there seems to be a dis- 
crepancy between observations and theoretical expectations, 
with observed lensing clusters being more centrally concen- 
trated than predicted by CDM simulations. This might also 
be the reason for some unexpectedly large observed Ein- 
stein radii (Zitrin et al. (2011)). Strong lensing bias might 
explain these differences (Comerford & Natarajan (2007); 
Oguri ct c. ( '■'); Meneghetti ot al (2011)) due to projec- 
tion effects for triaxial halos, baryons, or of foreground or 
background objects. 
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At present, strong gravitational lensing is the only avail- 
able tool to detect substructures within the lensing mass 
distribution beyond the Local Universe and independently 
of the luminous content. Faint and even dark substructures 
can be detected in lens galaxies either by observing multi- 
ply imaged lensed quasars with flux ratio anomalies (Mao 
& Schneider (1998); Dalai & Kochanek (2002); Metcalf & 
Zhao (2002); Kochanek & Dalai (2004); Bradac et al. (2002); 
Fadely Keeton (2012).) or via the gravitational imaging 
technique on Einstein rings and multiply imaged arcs (Veg- 
etti & Koopmans (zuUy); vcgctti et ai. (2UiU, 2Ui^)) While 
most of the observational and theoretical work that has been 
done to date is at the scale of galaxies, the current search for 
mass substructure can be extended to galaxy clusters using 
gravitationally lensed giant arcs. The large magnification of 
these arcs makes them sensitive to substructure masses as 
small as ^ 10^ — 10^ M© that can be detected using the 
gravitational imaging technique. 

Constraining the fraction of mass substructure at the 
scale of galaxy clusters is also important for understanding 
the properties of high redshift galaxies. Many of the proper- 
ties (e.g. stellar mass and star formation rate) of high red- 
shift galaxies detected using clusters as cosmic telescopes, 
can be measured accurately only when the source magnifi- 
cation is known. However, mass substructure introduces sig- 
nificant fluctuations in the source magnification and needs 
to be properly accounted for. In principle, numerical sim- 
ulations could be used to quantify the effect of undetected 
substructure on the source magnification. In practice, how- 
ever, the resolution of the simulation could potentially lead 
to an underestimate of the substructure fraction in the inner 
regions, while particle noise could mimic the effect of sub- 
structure and introduce spurious effects. In general, from 
a numerical point of view, when comparing theoretical pre- 
diction from numerical simulations and results from observa- 
tions, one should carefully quantify the effect of the limited 
resolution of the simulation and in particular, the effect of 
the particle noise. 

In this paper, we use the highest available resolution 
cluster simulations, the Phoenix simulations by Gao et al. 
(2012) to simulate gravitational lens clusters. Our main goal 
is to investigate the properties of the particle noise with 
a focus on multiple lensing properties and gravitationally 
lensed images. The first part of the paper is focused primar- 
ily on the quantification of the particle noise for numerical 
simulations of gravitational lensing, while the second part 
is focused on comparing the effects of particle noise with 
those of physical substructures. In particular, the paper is 
organised as follows. In Section 2 we briefly summarize the 
details of the Phoenix simulations. In Section 3 we will intro- 
duce the numerical lensing methods. We then quantify the 
effects of the particle noise on results from N-body simula- 
tions in Section 4, considering major lensing quantities such 
as the surface mass density, the deflection angles, the shear, 
the magnification, the critical lines and the lensed images. 
Section 5 investigates the scaling of the noise with number 
of particles and smoothing scale. In Section 6, we compare 
the effects of noise with those of physical substructures on 
each of the lensing properties. This comparison enables us to 
calculate a resolution limit for the smallest detectable sub- 
structures in a N-body simulation of gravitational lensing as 
a function of the particle number in the simulation. 



Table 1. Properties of the cluster E of the Phoenix simulations. 

Level mp[M0] A^200 e[kpc] 

4 1.39 X 10^ 5.8 X 10^ 3.84 
2 6.06 X 10^ 1.3 X IQS 0.44 



2 NUMERICAL SIMULATIONS 

In this paper we use the Phoenix simulations by Gao et al. 
( ) which are resimulations of 9 clusters from the Mille- 
nium Simulation ( ^rinpol ot al ("^OO.^)) . The simulations 
span masses from 7.5 X 10^^ M© for cluster C to 3.3 x 10^^ M© 
for cluster I within virial radii of 1.9 Mpc and 3 Mpc re- 
spectively. All Phoenix clusters are simulated at two reso- 
lution levels, Level-2 and Level-4, except for one cluster for 
which two additional resolution levels, Level-3 and Level- 1, 
are available. The 3D profiles of the simulated N-body clus- 
ters from the Phoenix simulations can be reasonably well fit 
by two components, a smooth cluster component and addi- 
tional small-scale substructures. The smooth component is 
triaxial and can be fit by either a NEW or an Einasto pro- 
file, although an Einasto profile seems to fit slightly better 
in most cases, with an average Einasto shape parameter of 
< a >~ 0.175. Eor a detailed description of the simulations 
and the parameters of all clusters we refer to Gao et al. 

( )■ 

In this paper we focus, as an example, on the cluster 
E, its properties are listed in Table 1. Cluster E is sim- 
ulated at two different resolutions and has a virial mass 
of M200 8.1 X lO^^Mo inside ^200 ^ 1-9 Mpc. Gao 
et ai. (20i^) quantify the performance of the NEW and the 
Einasto fit with a figure of merit function Q^, defined as the 
squared logarithmic deviation (Inp — lnp™°^^^)^ averaged 
over logarithmic radial bins. Eor example for the E halo an 
Einasto profile with Q — 0.067 fits considerably better than 
a NEW profile with Q = 0.135. Gravitationally bound sub- 
structures (subhalos) in the simulation with more than ^ 20 
particles are identified with the SUBEIND algorithm (for a 
recent comparison of different subhalo finders see jus 
et al. (2012)). Gao et al. (2012) describe the properties of 
the subhalo population in detail. In the mass range from 
10~^ < Msub/M2oo < 10~^ the number of substructures in- 
creases with decreasing subhalo mass as dN/dM oc M~^'^^ . 
The spatial subhalo distribution shows a distinct central core 
and the subhalos are found mainly at larger distances from 
the centre of the cluster. 

Eor the lensing simulations throughout this paper we 
use the same cosmological parameters as in the simulations, 

= 0.25, Qa = 0.75 and h = 0.73. The N-body simulation 
stores outputs for redshifts uniformly spaced in log a, where 
a is the scale factor. Eor the analysis presented in this paper, 
we consider the snapshot at redshift z — 0.32 and we place 
our simulated sources at redshift z = 2.0. These values are 
chosen such that the cluster is a relatively efficient lens, and 
the lensing configuration is comparable to observed strong 
lensing clusters, for example, the median cluster redshift for 
the CLASH survey is 2; = 0.4 (rostman ci, eta. (^ui.)), while 
the mean source redshift distribution of the SLOAN giant 
arcs survey peaks at 2; = 1.821 (Bayliss et al. (2011)). 
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3 LENSING THEORY 

In this section, we describe how we process the simulation 
data to get a high resolution simulation of the gravitational 
light deflection by the N-body cluster. 

During the N-body simulation, snapshots of the positions, 
velocities and other parameters of the particles were stored 
for multiple redshifts. We use these stored outputs of the N- 
body simulation for an accurate numerical simulation of the 
gravitational lensing effect. In particular, we use the 3D po- 
sitions of the particles of the N-body simulation to calculate 
a smoothed 2D mass density distribution. We then use this 
smoothed particle distribution to calculate the main lensing 
properties such as the deflection angles a, the magnification 
/i, the shear 7 and the lensed images. 

We adopt the smoothing algorithm commonly used in 
smooth particle hydrodynamics (SPH) simulations (Mon- 
aghan (1992)) to get a smooth density distribution S (x) 
from the A^part particles of the N-body simulation. We as- 
sign a different smoothing size Ip to each particle. This par- 
ticle size is variable and adapted to the local mass density 
in 3D. More specifically, this smoothing size is chosen to be 
equal to the distance to the A^ngb-th neighbour in 3D. By 
doing so, we use the full three dimensional density infor- 
mation to calculate the smoothing lengths of the particles. 
This is computationally more expensive than the equivalent 
2D adaptive smoothing, but it results in a density map with 
an enhanced contrast (Li et al. (2006)). The exact choice 
for A^ngb depends on the form of the kernel. Increasing the 
number of particles in the kernel will increase the smoothing 
and will therefore reduce artificial particle noise, but it will 
also smooth physical substructures. In this paper, we use 
a fourth order polynomial and the distance to the nearest 
64 neighbours, for more details see Sec. 4.2. In Sec. 5, we 
will investigate the dependence of the noise on the number 
of neighbours. Once the smoothing length for each particle 
is calculated, we use two different methods to simulate the 
gravitational lensing by this smoothed N-body mass distri- 
bution. 

The first method is based on a discretization on high 
resolution grids. In order to calculate the surface mass den- 
sity S (cc), all particles are projected onto a 2D grid using a 
smoothing kernel, of smoothing length chosen as described 
above. The surface mass density is then scaled by the critical 
density Ecrit = Ds / {^nGDdDds) to get the convergence 
HI (cc) = S (cc) /Scrit whcrc Ds, Dd and Dds are the angular 
diameter distances from the observer to the source, to the 
lens and from lens to source respectively. Finally, the gravi- 
tational lensing potential ^ (cc) (Bartelmann (2003)) is cal- 
culated by convolving the scaled surface mass density k (cc) 
with a logarithmic kernel, ^ (cc) — ^ j d^x' [x') In |cc — cc^|, 
using a Fast Fourier Transformation (FFT). Once the po- 
tential has been obtained, the lensing quantities a, 7 and /i 
can be calculated by using derivatives of the 1st, the 2nd 
or combinations of the 2nd order, respectively. The deriva- 
tives can be done numerically via finite differencing, or in 
Fourier space directly. For the highest resolution simulations 
we project the particles onto two grids of 16384^. Because 
of zero padding, the size of both grids is increased by a fac- 
tor of two by the FFT. The coarse outer grid covers a large 
area of (40 • 230^^)^ around the cluster and includes all ex- 
ternal shear components. The fine inner grid covers an area 



of (230^^)^. This corresponds to a maximum resolution of 
0.014^Ypix in the inner parts of the halo, which is a factor 
of 5 smaller than the softening length of the original sim- 
ulations at Level- 2. Both grids are centred on the cluster. 
In order to reduce the computational cost for the bootstrap 
resamplings described in Sec. 4, we reduce the grid size to 
2048^ 

The second method for calculating the lensing quan- 
tities is the smooth particle lensing method (SPL) intro- 
duced by Aubert et al. (2007). This method uses Gaussian 
shaped particles with a smoothing length Ip and takes ad- 
vantage of the linearity of the lensing quantities, X (r) = 
Ep^r^pW, whereXG{a , 7, /i, . . .}. As an example the 
deflection angle at any point on the lens plane is the sum 
of the contributions of the deflections of all particles of the 
N-body simulation. This second method is equivalent to con- 
volving the 3D particle distribution with a specific kernel for 
each of the lensing properties (Aubert et al. (2007)); for the 
deflection angles in the r direction, for example, the kernel 
reads as follows, 



exp(^) 



(1) 



where r is the distance of each particle p to the evaluated 
point. The direct computation of the convolution is very ex- 
pensive since it requires the evaluation of A^part x A^ray values. 
However, since the kernel depends on the distance as 1/r, 
there exist methods in order to reduce the computational 
load. The most popular solution makes use of a tree code to 
sort and group all particles at large distance and approxi- 
mate their contribution to any lensing property at a given 
point with a particle of bigger mass (Aubert et al. (2007)). 
In the rest of the paper, we will mostly use the FFT method, 
but we will also make use of the SPL formalism to quantify 
the particle noise analytically. 



4 PARTICLE NOISE 

The mass distribution of any N-body simulation is dis- 
cretized by individual mass particles. This discrete numer- 
ical sampling of the underlying mass density distribution 
introduces shot noise, which will affect all the major gravi- 
tational lensing quantities. In particular, the noise level will 
depend on parameters such as the N-body particle number, 
the particle size and the smoothing algorithm that is used 
to process the N-body particles. In this paper we calculate 
the noise for Poisson distributed particles. We are therefore 
assuming that each particle is sampled from a Poisson distri- 
bution with mean and variance one and that this sampling 
is equivalent to the Poisson sampling of an underlying true 
mass distribution. 

The theoretical basis for the expectation value and the 
covariance analysis for the interpolation of irregular sam- 
pled measurements to a smooth map is covered in a series 
of three analytic papers by Lombardi & Schneider (2001, 
2002, 2003). Focusing on gravitational lensing, Amara et al. 
( ) study the effect of particle noise for an analytic sin- 
gular isothermal ellipsoidal mass model and for a N-body 
cluster. They add artificial 2D Poisson noise to their sim- 
ulated surface mass density and then visually compare the 
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noise in the critical curve and caustic and the magnifica- 
tion of lensed images for several constant smoothing kernel 
sizes. Bradac et al. (2004) estimate the mean noise on the 
projected mass density within the critical lines by using 10 
bootstrap resamplings (compare Sec. 4.1) of the particles 
of their numerically simulated cluster. They, however, use 
a smoothing algorithm based on a Delaunay tessellation of 
the N-body particles for their lensing simulation. For their 
analysis of the cusp relation, a relation between the magni- 
fication of multiply lensed images close to a cusp (Schneider 
vvcisto (x^^^)) they also compare the particle noise for the 
Delaunay tessellation smoothing with a more simple Gaus- 
sian kernel smoothing algorithm with fixed smoothing size. 

Li et al. (20^ ) suggest adapting the size of the smooth- 
ing kernel in order to improve the contrast of the lensing sim- 
ulation. They estimate the noise of their 3D density adaptive 
smoothing algorithm for a uniform density field, for simu- 
lated isothermal ellipsoids and for a numerically simulated 
cluster and compare the results with the results from the 
simpler 2D adaptive smoothing. They simulate the particle 
noise by randomly populating fitted elliptical contours of 
the smooth halo with the same number of particles as the 
original simulation. They conclude that most of the wiggles 
in the critical line are not due to substructures in the sim- 
ulation. Their simulations of the cusp-caustic relation for 
the N-body cluster and the Monte-Carlo re-sampled cluster 
show that the particle noise produces many high-order sin- 
gularities. They also show that a more elliptical projected 
density is more sensitive to high-order singularities of the 
caustic. This is particularly important for simulating the 
anomalous flux-ratio problem with N-body simulations. 

Although several authors have addressed the problem of 
the discreteness noise of a N-body simulation and have es- 
timated the effect on some lensing properties, to our knowl- 
edge, to date no one has systematically investigated the im- 
plications for the simulation of gravitational lensing for very 
high resolution simulations. 

In this paper we use the smoothing algorithm from 
et al. (200ij) and analytically and numerically investigate 
the magnitude of the effects of particle noise on the lensing 
properties. We extend the investigations of Li et al. (2006) on 
the particle noise with a focus on different lensing properties 
and multiply lensed images and compare the noise to the 
simulated substructures for a state-of-the-art N-body galaxy 
cluster simulation. In this section we first present the two 
methods that we use to quantify the particle noise in the 
simulation; we then discuss in more detail how the particle 
noise affects the individual lensing properties such as the 
surface mass density, the magnification, the deflection angles 
and the lensed images. 

4.1 Bootstrap 

In order to simulate the noise numerically, we create B boot- 
strapped resamplings of the particles from the N-body sim- 
ulation. Each bootstrapped resampling is created by ran- 
domly choosing A^part particles with replacement from the 
A^part particles of the simulation. In this way, some particles 
will be included more than once, while others will not be in- 
cluded at all. We then project all particles on a grid and cal- 
culate the lensing potential, the deflection angles, the shear 
and the inverse of the magnification for each of the B re- 



samplings. Because this calculation is very time consuming, 
we restrict our calculation to B — 100. This is sufficient to 
calculate statistically converged results for the lensing prop- 
erties that we consider in Sec. 4.2 to 4.5. 

To test whether the bootstrapping is a good estima- 
tor of the inherent particle noise, we test the bootstrapping 
method as follows. We randomly chose 0.01 A^part particles 
from the high resolution Level- 2 simulation with replace- 
ment. This factor should be large enough to sample the 
total mass distribution randomly, independently of corre- 
lated phases, caused by the dynamical evolution during the 
N-body simulation. We repeat this 50 times, and then cal- 
culate the variance of those 50 randomly sub-sampled clus- 
ters. To test the validity of the bootstrap method, we then 
randomly choose one additional subset of N = 0.01 A^part 
particles from the original A^part particles with replacement 
and multiply bootstrap it. We, therefore, randomly chose 
50 times A" of those last A" particles with replacement. The 
two noise estimates for the subset of 0.01 A^part are identical. 
Therefore we can use the bootstrap method to estimate the 
particle noise of the original simulation. 

For a high resolution lensing simulation, calculating the 
noise via the bootstrapping technique can be computation- 
ally very expensive. In particular, most of the computational 
effort is spent in projecting the simulation particles onto the 
two-dimensional grids. In the following, we present, there- 
fore, a new way to calculate the noise analytically. 

We will make use of the SPL method ( Aubert et al. 
(900'^)) to derive an analytic expression for the bootstrap 
noise. The SPL method evaluates the contribution of each 
particle separately, and it is, therefore, ideal to evaluate the 
noise caused by the discreteness of the N-body simulation. 
We assume uncorrelated and Poisson noise for the 3D parti- 
cle distribution, that is, we assume for every particle a Pois- 
son probability density distribution with a mean and vari- 
ance of one of being included in the resampling. For a par- 
ticle convolved with a kernel, the Poisson noise is smoothed 
out over the size of the kernel. The total variance at any 
point X on the lens plane is therefore a sum over all the 
uncorrelated variances of the individual N-body particles, 

iVpart 

a^(x)= W^{\x-x',lk), (2) 

where Y G {a^, a,7} and VKy(|cc — x'\,lp) is the appropri- 
ate SPL kernel listed in the following sections. This method 
calculates the noise directly without the need of numerical 
bootstrapping. However, evaluating Eq. (2) numerically for 
quantities like for example the deflection angles, is labori- 
ous. The kernel for the deflection angles is not compact and 
therefore for any point on the lens plane x one has to sum 
over the long-range contributions for all N-body particles. In 
order to effectively evaluate Eq. (2), one has to use some kind 
of approximation. Two possibilities are, either group distant 
particles and approximate the effect by a bigger particle in 
a tree code, or to use a two-dimensional approximation (see 
Sec. 5.2). 

4.2 Surface Mass Density 

The most basic quantity that can be derived from the dis- 
crete N-body particle distribution is the smoothed surface 
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mass density. In the grid-based approach we derive all other 
lensing properties via FFT from the projected particle dis- 
tribution, see Sec. 3. In this section, we describe therefore 
in detail, the smoothing algorithm and the noise properties 
of the projected mass density distribution. Simulations of 
gravitational lensing for a long time used fixed Gaussian ker- 
nels to smooth the particles of the N-body simulation (e.g. 
Bartelmann et al. (1998); Bradac et al. (2002); Maccio et al. 
(2006)), or adapted the smoothing size if too few particles 
fell into the kernel (e.g. Li et al. (2005)). A more sophis- 
ticated method was then proposed by Bradac et ai. (2004) 
who used a Delaunay tessellation (Schaap & van de Wey- 
gaert (2000)) to obtain a smooth density distribution from 
the N-body particles of a numerical simulation. T i et al. 
(2006) suggested to adapt the size of the smoothing to the 
3D density distribution in order to improve the contrast of 
the lensing simulation. In this paper we will use the method 
from Li et al. (2006), and use the fully 3 dimensional particle 
distribution to determine the size of the smoothing kernel. 

The 2D surface mass density S (cc) is calculated from 
the A^part particles of the N-body simulation. Each particle is 
projected with a normalized kernel, f d^x'Wi{\x—Xi\, U) — 1 
in order to obtain a smooth density distribution. 



(3) 



The size of the kernel, is different for each particle i and is 
calculated from the distance to the nearest A^ngb neighbours 
in 3D. 

There are different forms of smoothing kernels, the 
most widely used 2D functions are Gaussian, Wi{r) — 

1/ {2iil{^ exp r^/ ^2/~f^j, a simple polynomial, Wi{r) — 

3/ (vr/f) (1 - r'^/li) for r < h, or the more complicated 
polynomial (Li et al (2006); Springe! (2005)), 




< r < |, 

1 <r </^, 
otherwise. 



(4) 

For a Gaussian smoothing width U = y^l03/1120 /^, the ef- 
fective area covered by a particle is equivalent to a particle 
smoothed by the two other kernels. The simulation of gravi- 
tational lensing does not strongly depend on the form of the 
kernel. Since the projection of the N-body particles on the 
lens plane is the most time consuming step, we choose the 
second kernel, which is the fastest to evaluate numerically. 

The convergence, the surface mass density in units of 
the critical density, is shown in the top panel of Fig. 1 for the 
Phoenix cluster E at Level-2 resolution. The figure shows the 
central part of the halo, the side length is about three times 
the size of the critical lines. The extent of the critical line 
is comparable to the k = 0.7 contours, the figure therefore 
covers the whole region where multiply lensed images occur. 
As an example, a second contour at At ^ 0.5 is over plotted. 
Both contours are not smooth and show several irregulari- 
ties. These 'wiggles' arise from two different effects. The first 
effect is due to the presence of physical mass substructure 
in the simulation. For example, in the central (92^^)^ region, 
shown in Fig. 1, there are 2597 subhalos identified by SUB- 
FIND with ATpart > 20, of which 2131 have a Msub < lO^M©, 
58 have a Msub > lO^^M© and 408 have an intermediate 




-40" -20" 



Figure 1. Top panel: Projected convergence of the N-body parti- 
cle distribution at Level-2 resolution for cluster E. The smoothing 
is chosen to be equal to the 64th neighbour in 3D. Convergence 
K for the central (0.55 Mpc)-^, the critical lines are approximately 
the same size as the contours dX k, = 0.7. The contours are ir- 
regular because of substructures and the effect of discreteness 
noise. Bottom panel: Particle noise of the convergence c (cc) in 
percent of the projected scaled surface mass density from 100 
bootstrapped resamplings of the particle distribution 



mass. The second effect, which is clearly visible in the top 
panel of Fig. 1, is related instead to the presence of parti- 
cle noise. Even for this very high-resolution simulation both 
effects are significant and comparable in the inner regions, 
r < 20'' , of the halo. In the following, we will quantify the 
contribution to the surface mass density fluctuations due to 
particle noise. 

By substituting one of the three smoothing kernels, for 
example Eq. (4), in the expression for the analytic variance, 
Eq. (2), we get an analytic expression for the particle noise 
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Figure 2. Critical lines for different particle noise realisations 
of the cluster. Background in grey scale is the particle noise of 
the inverse of the magnification 3(j^-i The noise on the 

inverse of the magnification is cut-off at 3a^-i ~ (white 
band following the critical curves), this indicates the noise on the 
critical line. 



on the surface mass density, 



(5) 



The bottom panel of Fig. 1 shows the noise (x) on 
the surface mass density as a function of hi{x). There is 
a significant difference to a calculation with 2D adaptive 
smoothing. In the 2D case, all line-of-sight N-body particles 
for a point on the lens plane are assigned the same smooth- 
ing length /. For this smoothing we expect the variance of the 
Poisson noise on the surface mass density to be proportional 
to the number of line-of-sight particles at each point. Since 
the number of particles can be estimated from the conver- 
gence K via A^Ecrit/^ (cc) /'Trip, we would expect the bottom 
panel of Fig. 1 to be roughly constant. In contrast, here, we 
adapt the smoothing to the local density in 3D, and there- 
fore, the noise depends on the integrated particle density 
distribution along the line of sight and is more complicated 
than the simple 2D case. For the Level-2 resolution of the 
Phoenix simulations in Fig. 1 the mean noise on the surface 
mass density is about 1 — 2% and < 5% even in high density 
regions. 



4.3 Shear and Inverse of the Magnification 

In the last section we quantified the noise on the surface 
mass density. We will now study the noise on other lensing 
quantities, namely the shear and the inverse of the mag- 
nification, by numerically calculating the lensing potential. 
Since the magnification and both shear components are im- 
portant for weak lensing studies, we will also discuss the 
analytic noise calculation from the particles of a N-body 
simulation. The shear and the scaled surface mass density 
are related to second derivatives of the lensing potential by 
71 = (^,11 - ^,22)/2, 72 = ^,12 and k = (^,ii + ^,22)/2, 
where ^,12 = dxdy^\ from these the magnification can be 
calculated as /i = [(1 — k)^ — 7^] 

Using the SPL kernel for the shear from Aubert et al. 
(2007) we get the kernel for the variance of the two compo- 



nents of the of the shear 7 = (71,72), 



ic -y ^ 



(6) 

where G{rp,lp) = exp [-rV (2/p)] and = + . We 
use this kernel to calculate the variance of the two shear 
components, cr^ = "^il^^ ^cr^ {xi — x^U). Due to the term 
C^, the support of the kernel (6) is more compact than the 
Gaussian smoothing kernel. Therefore Eq. (6) can be easily 
evaluated by direct numerical summation. This provides a 
fast and easy way to calculate the particle noise on the two 
shear components. 

Taylor expanding the inverse of the magnification allows 
us to derive the noise on the inverse of the magnification due 
to noise on 7 and 



4(1 



\2 2 



. A 2 2 
+ 471 CT^, 



I /I 2 2 
+ 472Cr^2 



(7) 



As an example. Figure 2 shows the 3a particle noise on 
the inverse of the magnification Sa^-i as a grey scale 

background. In order to enhance the effect, the noise level is 
set to white for all pixels where 3cr^-i > close to the 

critical line. This indicates the width of the 'wiggles' of the 
critical line due to discreteness noise in the region close to 
the critical line where 0. The critical lines for four 

random particle noise realisations of the numerical cluster 
are over plotted with coloured lines. The blue line is for the 
original N-body cluster E from the Level-2 Phoenix simu- 
lations. The discreteness of the particles in the simulation 
is especially important in the strongly magnified and highly 
nonlinear regime of the critical lines. Figure 2 demonstrates 
the magnitude of the effect that significantly shapes the crit- 
ical lines, even for this high-resolution simulation. Multiple 
imaging caused by strong lensing also takes place in this 
highly magnified region. It is therefore important to under- 
stand and quantify the noise on the deflection angles and the 
highly magnified images. We describe the properties of the 
particle noise on the deflection angles in the next section. 



4.4 Deflection Angles 

Understanding the noise on the deflection angles is essential 
to quantify and describe the noise on the lensed images. We 
calculate the noise on the deflection angles by projecting all 
particles with a modified kernel, 



2G'(rp,/p)+ 1 



(8) 



Here W^2 {x) is a vector for the noise on the deflection an- 
gles in the x- and y directions, c = mp/(7rEc), Vp is the 
distance of particle p to the point CC , Ip IS the smoothine: 

length and G(rp^lp) — exp(-2y^). The first two terms in 

A can be easily calculated by direct numerical summation. 
This is necessary since the smoothing length Ip is different 
at each point. For the summation it is sufficient to consider 
points with Vp ^ 5/p, since at larger distances for both terms 
in A, the Gaussian kernel vanishes, G{rp,lp) 0. The sec- 
ond term, B, is easily evaluated by convolving the grid via 
FFT. 

The noise on the deflection angles in x direction is shown 
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Figure 3. Difference between two random realisations of the 
particle noise of the deflection angles in x direction, — (left 
panel). The marked region shows the size of a typical image from 
Sec. 4.5. Standard deviation of the deflection angles in the marked 
region calculated with Eq. (8) (right panel). For both panels the 
colour scale is in arcseconds. 

in the right panel of Fig. 3 for the Level-2 version of cluster 
E. The side length of the figure is comparable to the size of 
a strongly lensed image. Over the size of a typical giant arc, 
the effect of shot noise on the deflection angles is a smooth 
function. The left panel of Fig. 3 shows a realisation of the 
particle noise on the deflection angles in the x direction. The 
marked region corresponds to the size of a typical multiply 
lensed image. Due to the 1/r dependence of the deflection 
angles, the particle noise is strongly correlated on long scales 
in the left panel of Fig. 3. These long-scale correlations will 
have a big effect on the calculation and the understanding of 
the noise on the lensed images presented in the next section. 

4.5 Lensed Images of a Gaussian Source 

In order to understand the effect of particle noise on 
the strongly gravitationally lensed images of a background 
source, we simulate a gravitational lens system by placing a 
source galaxy at z = 2.0 behind cluster E at Level-2 resolu- 
tion, (4.4^^,4.1^^) away from the cluster centre in projection. 
The source has a Gaussian surface brightness distribution 
with a FWHM size of 0.51^^ and a peak surface brightness 
of 100. The resulting gravitationally lensed images are shown 
in the top left panel of Fig. 4. Multiple images of a back- 
ground source form in regions of high magnification, where 
the effect of the particle noise, due to the discreteness of 
the N-body simulations, is the largest. It is very important, 
therefore, to understand how the noise on the cluster sur- 
face mass density and deflection angles is propagated to the 
lensed images. In the following, we quantify this noise from 
the bootstrapped re-samplings of the N-body particles (com- 
pare with Sec. 4), as well as from analytical approximations. 

For each bootstrapped particle distribution i of the clus- 
ter lens, we calculate the deflection angles cXi and use them 
to lens an identical source surface brightness distribution s 
for each bootstrap resampling of the lens. The correspond- 
ing images, di, differ from each other because of the particle 
noise. 

The lower left panel of Fig. 4 shows the standard devi- 
ation of 100 lensed images of the same background source 
calculated in the same way from 100 different bootstrapped 
cluster resamplings; the noise on the image brightness is as 
large as ^ 10% of the image surface brightness. The ori- 




Figure 4. Image of a Gaussian source (top left) using the N-body 
cluster, brightness difference of two bootstrapped realisations (top 
right), standard deviation of 100 resamplings (bottom left) and 
analytic approximation using the noise of the deflection angles 
and a linearization of the lensing equation (bottom right). 

gin of the image brightness fluctuations due to the particle 
noise and their distribution can be understood qualitatively 
by considering the lens equation y = x — cx (x). According 
to this equation, each point x on the lens plane is deflected 
by a deflection angles a to a point y on the source plane. 
As shown in Sec. 4.4, the noise associated with the discrete 
particle distribution is responsible for local fluctuations on 
the gravitational potential (x) of the lensing cluster and 
therefore on the deflection angles a = V'ip{x). The noise 
on the lensed images is then just the effect of the particle 
noise on the cluster surface mass density propagated to the 
images via the deflection angles. From the left bottom panel 
of Fig. 4, it is clear that a large number of realisations is 
needed in order to reduce the statistical error of the boot- 
strap noise estimation. More quantitatively, we can extend 
the analytic noise analysis of the deflection angles in Sec. 4.4 
to the noise on the simulated images, by linearizing the lens 
equation and approximating the difference in the surface 
brightness between two images 6d by (ivoopmans (2005)), 

6d ^ —{6cxx dxS -\- 6cxy dys). (9) 

For a given source surface brightness distribution, fluctua- 
tions in the image surface brightness are a linear combina- 
tion of the fluctuations in the deflection angles. The bottom 
right panel of Fig. 4 shows this approximation for the noise 
on the images using the noise on the deflection angles. If 
we compare the two lower panels of Fig. 4, we see that this 
linearization overestimates the noise where the noise is high. 
From the above equation, we also expect the noise to follow 
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Figure 5. Caustics of two different particle noise realisations. 
Due to long-scale correlations of the deflection angles (cf. Fig. 3) 
and the flnite size of the image, random particle noise will result 
in a shift of the caustic. 



the distribution of the gradient of the source surface bright- 
ness. In particular, small changes in the deflection angles 
can be strongly amplified if the gradient of the source sur- 
face brightness distribution is large. In Sec. 4.6, we discuss 
the effect of particle noise for a source surface brightness 
distribution with varying level of smoothness. 

The top right panel of Fig. 4 shows the difference in 
brightness do — di for each pixel of the image plane, as an 
example, for two particular cluster re-samplings i = and 
i — 1. The observed surface brightness difference can be 
attributed mainly to two effects: A global shift to the left 
by < 0.1^^ of di relative to do, and smaller scale differences. 
We start by discussing the origin of this global shift. We 
will see in the following that this shift is equivalent to an 
unobservable shift of the caustic relative to the source. We 
will therefore present a method to separate and subtract 
the contributions of this shift to the noise on the image 
brightness later on in this section. 



4.5.1 Global shift of the lensed images 

Although, small scale fluctuations of the lensed images can 
be significant, an overall shift in the image position also 
contributes significantly to the difference between two boot- 
strapped images. In the following, we shall discuss the dif- 
ferences in the deflection angles and the caustic structure 
between two re-samplings in more detail. This will help us 
to eliminate most of the brightness difference in the top right 
panel of Fig. 4 by a relative source - caustic shift, and will 
provide a more physical and observationally motivated mea- 
sure for the noise on numerically lensed images. 

In the left panel of Fig. 3, the difference between two de- 
flection angle maps calculated from two bootstrapped clus- 
ters is shown. The size of a typical image, as for example 
in the top left panel in Fig. 4 is marked as a rectangle. Be- 
cause of the finite size of the area covered by an image on 
the lens plane and because of the long scale noise correla- 
tions on a, the mean of the noise realisation of the deflection 
angles within the area covered by the image is positive. This 



positive deviation of the deflection angles, also visible as a 
global shift of the caustic structure (see Fig. 5), will then 
result in a mean shift of the lensed images, or equivalent ly 
of the source on the source plane. Therefore, while a mean 
constant additional deflection angle does not change any of 
the physical parameters of the lens, it changes the relative 
position of the caustic and the (arbitrarily fixed) source. 

4.5.2 Small scale surface brightness fluctuations 

Any shift of the caustic position can be compensated by 
shifting the source position accordingly. Therefore, a sim- 
ple shift in the source position eliminates most of the image 
brightness difference in the top right panel of Fig. 4. In the 
previous section we therefore made an assumption. We used 
an identical source for the shifted and the not-shifted caus- 
tic structure and we therefore compared two non-equivalent 
images. 

Following a similar chain of argument, it is possible that 
the re-sampled cluster deflection angles will differ in higher 
order derivatives of the deflection angles as well, such as 
magnification, shear or flexion or more likely, a combination 
of all of those. Since none of the sources intrinsic properties 
such as position, size or morphology are known, we will as- 
sume in the following as little as possible about the source. 
The method described below is closer to an observational 
point of view and only constrains the source by its regular- 
ity. We therefore rephrase the problem of the comparison of 
two equivalent images as follows. Two deflection angle maps 
are given and the first of the two images is fixed as reference 
image. What is the most similar image we can find using the 
second deflection angle map, if we are only allowed to use a 
relatively regular source? 

This question can be more easily answered within a 
Bayesian formulation of gravitational lensing. We try to find 
the closest image di, to a reference image do under some 
regularity conditions. For a given, fixed smoothness of the 
source it is in general not possible to perfectly recreate the 
input image do exactly using the deflection angles ai and 
a smooth source Si. It is, however, possible to reconstruct 
a close image di . We will try to find the best reconstructed 
image di, which is as close as possible to the input image 
do, and at the same time keeping its source Si regular. Our 
source reconstruction is pixelized on a non-regular triangu- 
lation and the smoothness is only constrained by a curvature 
regularisation (Vegetti & Koopmans (2009)). For a given in- 
put image do we maximize the Bayesian posterior 



\ogP{si\do,OLi,\Hs) ^ X + K\\HsSi\\i (10) 

where is between both images, Hg is the form and As 
the strength of the regularisation of the source. We assume 
a quadratic Gaussian prior for the regularisation which is 
centred at s = 0. The regularisation strength As is found 
self consistently from the data themselves by maximising 
the posterior 



P(As|do,ai,Jfs) 
oc j dsi P(do|si 



^ P(do|As,ai,Jf.)P(As) 
P(do|ai,ifs) 

As,ai,Jfs)P(5i|ai,Jfs). 



(11) 



Here we assume a prior P{Xs) which is flat in log As. This 
method uses minimal assumptions about the source to find 
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Figure 6. Difference in the image brightness using a Bayesian 
source reconstruction (top left panel), colorscale identical to Fig. 4 
top right panel. Difference between original source and recon- 
structed source shows that the biggest effect is a shift in position 
(top right panel). Standard deviation of 100 bootstrapped real- 
isations of the N-body cluster, for each cluster the best image 
matching an input image is found (bottom left panel). Standard 
deviation approximation with analytic formula using ctq, and a 
linearization of the lensing equation (bottom right panel). 



the best matching image for an input image. The pixelized 
brightness difference do (x) — di (x) of the reconstruction 
is shown in the top left panel of Fig. 6 for the same ex- 
ample as in the top right panel of Fig. 4. This method 
automatically corrects for the shift in the deflection angles 
or in the caustic in Fig. 5 and compares images produced 
by equivalent sources. The surface brightness difference be- 
tween the two sources is shown in the top right panel of 
Fig. 6. As explained above, the reconstructed source is auto- 
matically shifted by the algorithm in order to find an equiv- 
alent source which now reproduces the reference image as 
well as possible. By applying this method to A^b = 100 
bootstrapped cluster resamplings we simulate the parti- 
cle noise on the pixelized image brightness distribution, 
al (x) = - 1) E^^^ (do - dif. The noise is shown 

in the bottom left panel of Fig. 6. By allowing the algorithm 
to automatically adapt the source, the noise on the images is 
reduced by a factor of two. In other words, half of the noise 
on the image brightness in the lower left panel of Fig. 4 can 
be attributed to a simple relative shift of the source to the 
caustic between different bootstrap realisations. The sharp 
cuts at the top and the bottom of the arc in the lower left 
panel of Fig. 6 are due to the caustic structure. The single 
imaged regions of the source are less noisy by a factor of 
^ 100 than the threefold imaged parts. The reconstruction 




Figure 7. Same as the bottom left panel of Fig. 4 but for a 
different source position and the source brightness distributions 
in the bottom panels. An increased structure of the source and 
therefore an increased source gradient increases the noise on the 
image brightness. 

algorithm for the source can more easily fit the single-imaged 
regions than the threefold imaged parts of the brightness dis- 
tribution. Therefore, the noise is not visible in these regions 
in Fig. 6. If we assume that the Bayesian source reconstruc- 
tion mainly corrects for the source shift, we also can approx- 
imate this noise analytically. Shifting the source keeps the 
gradient of the source Vs constant on the image plane. But 
we have to correct the deflection angles by the mean source 
shift, effectively calculating 5d ^ —{5ol — {5ol))V s where the 
mean (^a) is evaluated over the size of the images. The bot- 
tom right panel of Fig. 6 shows the corrected version of the 
linear noise approximation for the image brightness. Even 
though this is a very simplified and fast approximation to 
the noise, the result is very close to the accurately simulated 
noise. 



4.6 Lensed Images of a Realistic Source 

The noise on the surface brightness distribution of the 
lensed images only weakly depends on the image position, 
but it strongly depends on the source gradient. This can 
be understood from the analytical approximation in Eq. 9, 
since the noise on the deflection angles in the top right 
panel of Fig. 3 is very smooth in the region where strongly 
magnified images occur at x ~ Ibl5^^ As an example we 
show the different noise maps of the image brightness in 
the top panels of Fig. 7 for two more structured sources. 
The respective sources are shown in the lower panels. Both 
sources are located at (-4'', 1.38''), the FWHM size is 
the same as the Gaussian source in Fig. 4, 0.51'', and the 
maximum brightness is scaled to 100. The source on the left 
is a smooth version of the right source. The magnitude of 
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the noise in the top left panel is similar to the noise on the 
image for a Gaussian source. This can be seen by comparing 
Fig. 7 to the lower left panel of Fig. 4. The distribution of 
the noise follows that of the source gradient on the lens 
plane. The more irregular source in the right panel also 
shows a significantly increased noise. This dependence is 
easily understood from the linearization in Eq. 9. 



5 NOISE SCALING 

In the previous sections, we have calculated the particle 
noise due to the discreteness of the N-body simulation on 
different lensing properties. The results were calculated for a 
simulation with A^part particles and a 3D adaptive smoothing 
algorithm that adapts the smoothing length for each parti- 
cle to the distance of its A^ngb-th neighbour in 3D. In this 
section, we will derive how the particle noise in the simula- 
tion of gravitational lensing changes as a function of parti- 
cle number A^part and the number of smoothing neighbours 
A^ngb- Increasing the number of particles inside the smooth- 
ing kernel, A^ngb, will result in a more smoothed mass dis- 
tribution and therefore reduced noise on all of the lensing 
properties. It will, however, also smooth out the physical 
substructures that were part of the simulation. Increasing 
the particle number while keeping the number of neighbours 
constant will increase the mass resolution of the simulation 
and therefore also reduce the noise, at the cost, however, of 
the computational load. 

In principle, we could numerically simulate a change in 
A^part and A^ngb by resampling the mass density and then re- 
calculate the smoothing length for each particle and all of 
the lensing properties for each of the re-projections of the 3D 
particle distributions. In practice, we develop a method that 
allows us to convert the 3D adaptive smoothing lengths into 
2D smoothing lengths and derive analytic expressions for the 
scaling of the noise in 2D, while preserving the information 
of the 3D particle density distribution. 



5.1 Smoothing in 2D 

Until now we calculated the noise on all of the lensing 
properties from the 3D particle distribution and we used 
a smoothing length that was adaptive with the 3D particle 
density distribution. We will now introduce an approxima- 
tion that allows us to keep most of the 3D information, while 
performing the noise calculations in 2D. The 2D version ob- 
viously increases the speed of the calculations substantially, 
and it also allows us to derive the scaling properties of the 
noise on the surface mass density and the deflection an- 
gles with the particle number and the number of smoothing 
neighbours. 

By comparing Eqs. (5) and (3), we find that on average, 
on scales larger than the smoothing length k of each particle, 
the contribution of each particle i to the variance of the 
surface mass density, cr|, and to the surface mass density, 
H, differs by a factor of 



fi oc 



9 



(12) 



of the smoothing kernel. Here, we used the kernel W oc 
(1 — r^//^)^ for r < /, for details and other kernels see 
also Appendix A. Eq. (12) states that smaller particles from 
high-density regions contribute proportionally with a fac- 
tor of l/l^ to the variance of the surface mass density. We 
can use this information to derive a 2D approximation to 
the noise on the surface mass density. To this end, we de- 
fine an effective smoothing length, /eff (x) in 2D, which is a 
integrated average of all line of sight N-body particles. 



eff (x) 



1 1 



-1/2 



(13) 



This allows us to assign a single smoothing length to a par- 
ticular 2D position on the lens plane, and simultaneously 
taking into account the 3D density distribution. Note that, 
in the limit of an extremely high-resolution grid on the lens 
plane and a finite N-body particle number, not all positions 
will be occupied, and /eff = /• With this simplification we are 
able to approximate the variance of the surface mass density 
of Eq. (5), in terms of the 2D surface mass density S, 

9 m 



(x) 



5.2 Noise in 2D 



Stt ll^ (x) A 



(14) 



In this section, we introduce some useful expressions for the 
noise on the lensing properties using the 2D smoothing ap- 
proximation derived above. This will simplify the derivation 
in the next section and will allow us to calculate the noise 
on a high-resolution grid in a fast and easy way, by elimi- 
nating the need for a tree-based evaluation of the long-range 
terms of the deflection angles (e.g. the 1/r^ dependence of 
the variance). 

We start by considering the definition of the effective 
smoothing length for each point on the lens plane /eff (x) 
as presented in Eq. (13). This smoothing length defines the 
correlation length of the noise for each point in 2D. All we 
need for the first order covariance matrix is the amplitude 
of the uncorrelated noise as at each point on the lens plane, 
which can be calculated from the correlated noise in Eq. (5). 
In other words, the correlation introduced by smoothing the 
original N-body particles with a smoothing kernel of size /eff 
has decreased the uncorrelated noise to the expression given 
in Eq. (5). By deconvolving each point on the lens plane 
with a smoothing kernel of size /eff, the correlated noise, cr|, 
is increased again by 



(x) 



(15) 



as demonstrated in the Appendix A. In order to obtain the 
uncorrelated noise amplitude, as, we therefore have to in- 
crease the noise from Eq. (5). 

These considerations allows us to simulate the effect of 
the particle noise on the lensing properties with a 2D convo- 
lution. For the variance we obtain the following expression. 



ay (x) = J d?xa^{x' 
and for a single particle noise realisation 



(16) 



The constant factor 9/(57r) depends weakly on the form 



Ay(x) = yd^xi?(x ,as)W^Y(|x-x |,/eff), (17) 
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where Y G {k^cx,^}, Wy{\x — cc^|,/eff) is the appropriate 
kernel and R{x\aj:) is a random number drawn from a 
Gaussian distribution with standard deviation a^. Eq. (17) 
can be understood as placing uncorrelated particles of size 
/eff with random mass i?(cc^as) at the 2D positions cc^ 

5.3 Scaling with A^part and A^ngb 

With this effective 2D formulation of the adaptive 3D 
smoothing we can study the scaling of the noise. An in- 
creased resolution of a N-body simulation samples the den- 
sity with more particles of smaller masses. This increase 
in particle number also decreases the particle noise due to 
the finer discrete sampling. With an increased particle num- 
ber, the number density of particles in 3D is also increased. 
Therefore the size of the smoothing kernel, which here is 
the distance to the A^ngb-th neighbour, also changes. From 
Eq. (14) we calculate the scaling of the noise on the sur- 
face mass density between different resolutions k and j of a 
N-body simulation as 



ml_ ( r (x) 
\V {x) 



(18) 



If all particles have the same mass, then the fractional 
change in particle mass is equal to the inverse change in 
particle number, /m^ = A^p^rt /A^part- The change in the 
smoothing of the particles, /, can be estimated from the 
change in particle number and the change in the number of 
neighbours by 



(x) _ I 
P [x) ~ I 



ngb part 
ngb part 



(19) 



With these two scalings, the fractional change in the noise 
of the surface mass density in Eq. (18) simplifies to 



' part 
^ part 



1/3 



ngb 



2/3 



(20) 



This first result allows us to estimate the change in the vari- 
ance of the surface mass density with changing particle num- 
ber of the simulation and a changed smoothing length to 
convert the N-body particles into lensing properties. 
We can transform Eq. (14) to 
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c^l {^) 

E2(x) ^ Stt /2 (x) A S (x) ' 



(21) 



For the special case of a uniform density field in 3D with 
equally distributed particles, the mean projected 2D sur- 
face mass density is S = mpNpavtA/ L"^ , where L is the side 
length of the cube of equally distributed particles. For this 
configuration, the distance / to the A^ngb-th neighbour can 
be calculated from 



part 



1/3 



(22) 



Substituting these two relations into Eq. (21) we obtain in 
units of the critical density Ecrit , the following relation 



K (x) 



^ ^ ngb part 



(23) 



which is similar to what (lA et al. (200^ )) found by numerical 
fitting (their Eq. (4)). The exact value of the proportionality 
constant depends on the exact form of the kernel. 



6 COMPARISON OF THE PARTICLE NOISE 
WITH SUBSTRUCTURE 

In this second part of the paper we compare the small-scale 
fluctuations due to two competing effects, the particle noise 
and physical mass substructures (subhalos). In particular, 
we are interested in quantifying the limit at which the effect 
of physical mass substructure becomes comparable to that of 
particle noise, and the substructure is considered too 'small' 
to be 'visible' above the noise level. To this end, we need 
quantitative measures of the effect of a substructure on the 
lensing properties. Mainly, we need a metric to answer the 
following questions: which lensing property is best to look 
at in order to compare the substructure to the noise? Which 
property of the substructure is the best one to quantify how 
'small' a substructure is? When can a substructure be con- 
sidered 'visible' above the noise level? 

Substructures in our N-body simulation are identified 
as gravitationally bound objects consisting of more than 
A^min 20 particles. For each substructure we can measure, 
among other parameters, its mass, size, ellipticity, density 
profile and circular velocity. For a simplified substructure 
model such as a singular isothermal sphere (SIS), the lens 
strength b is proportional to the Einstein radius which is 
proportional to the one-dimensional velocity dispersion ay. 
The mass within R on the other hand is proportional to 
ayR. Therefore we choose the quotient of the mass of the 
substructure and its half-mass radius in units of the critical 
density, /3 = Mss/i^HMR/Scrit as a measure of the strength 
of each substructure. We will confirm in Sec. 6.3 that this 
measure derived from a simplified SIS model is a good pa- 
rameter to quantify the effect of the numerically simulated 
substructures. 

For the convenience of the reader we will also introduce 
a second x-axis on the top of Figs. 10 to 13. This second 
axis converts the lens strength on the bottom x-axis to a 
typical substructure mass, Msg. The conversion is a linear 
fit to all numerically simulated substructures as identified 
by SUBFIND within the Level-2 simulation of cluster E in 
log(Mss)-log(/3) space within —2 < log^^Q fiss < 0. This corre- 
sponds roughly to 0.03 < Mss/(1O^°M0) < 7. The linear fit, 
log(Mss/(lOi°M0//i)) = 1.21±o•OMoglo(As^) + 0.71=^°•°^ 
averages different substructure profiles, sizes and concentra- 
tions and therefore yields a typical substructure mass for a 
given substructure lens strength. All calculations, however, 
are performed in terms of the substructure lens strength /3 
in order to fully take into account the different 2D profiles 
of the simulated substructures. 



6.1 Substructure Surface Mass Density 

As a first, although very simple step, one might compare the 
projected scaled surface mass density of each substructure, 
Ki, to the amplitude of the noise fluctuations of the surface 
mass density due to the particle noise, a^. We approximate 
the noise by a Gaussian fluctuation field with amplitude 
and a correlation length equivalent to the 2D smoothing 
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Figure 8. Approximation of the lens strength /S of Gaussian noise 
with the same parameters as the noise on the surface mass density 
for the Level-2 version of cluster E. At the critical lines the lens 
strength of the 3cr noise is /3 ~ 0.14. 

length, /eff, at each point on the lens plane. We can then cal- 
culate the equivalent lens strength for the noise fluctuations. 
The size of the Gaussian is /eff = the total mass of a na 
noise substructure is therefore Trl^na-s for n G {1, 2, 3}. The 
half mass radius is i?HMR = /V2 In 2. Therefore the equiva- 
lent lens strength is /3 — Mss/{RuMR^crit) = 7r/crK/(2 In 2). 
Figure 8 shows the lens strength of the noise on the scaled 
surface mass density at each point on the lens plane. At the 
critical line, we therefore identify the substructures that will 
have a smaller influence than the 3a noise as those substruc- 
tures which are smaller than /3 ^ 0.14, this corresponds to 
an average subhalo mass smaller than Mss ~ 6.5 x lO^M©. 

By comparing the surface mass density, we are essen- 
tially comparing a single substructure with a random field 
of positive and negative substructures for the noise. Since 
the lensing equation is highly nonlinear, especially in the 
strong lensing regime, the noise also propagates nonlinear ly 
through the lensing properties as we will show in the follow- 
ing sections. We can therefore not assume that the limits 
derived here based on the surface mass density are the true 
resolution limits of the simulation, in the sense that every 
lensing property yields identical limits. 

6.2 Substructure Magnification 

As a second quantity we compare the effect of mass sub- 
structure and particle noise on the inverse of the magnifi- 
cation and on the critical lines. Both panels in Figure 
9 show the particle noise on the inverse of the magnifica- 
tion 3cr^-i/|/i~^ I as a grey scale background. The noise in- 
creases with the magnification and reaches its maximum at 
the critical lines. In order to indicate the width of the 3a 
'wiggles' in the critical line, values where the noise 3a ^-i 
exceeds close to the critical curves are set to white. In 

the top panel, two critical lines, ^0, are over plotted. 
The blue line is for the original N-body cluster E from the 
Level-2 Phoenix simulations. The red line is for the same 
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Figure 9. Top panel: Critical lines with (blue) and without (red) 
subhalos. Background in grey scale is the particle noise of the in- 
verse of the magnification 3a^-i /\iJ,~-^\. The noise on the inverse 
of the magnification is cut-off at 3cr^-i ~ Im""*^! (white band 
following the critical curves), this indicates the noise on the crit- 
ical line. Coloured circles indicate substructures with masses in 
1O-'^^M0. The inset shows an enlargement with the significant ef- 
fect of three (four) subhalos directly on the critical line marked 
by arrows, for details see text. Bottom panel: 11 substructures of 
decreasing lens strength artificially added on top of the N-body 
cluster mass distribution at — 15.6^^ —2.6". Substructures with a 
lens strength smaller than ~ 0.1^^ are within the 3a limits of the 
noise on /i"-*^ shown as a grey scale background. 



cluster but with all subhalos identified by SUBFIND re- 
moved. The physical subhalos of the simulation are marked 
as coloured circles. The colour indicates the subhalo lens 
strength /3 = Mss/i^HMR/ScHt- The two critical lines are al- 
most identical except for those few cases in which a subhalo 
lies directly on top of the critical line. The inset is show- 
ing an enlargement of three of these cases indicated by ar- 
rows. The masses of the three subhalos are 6.0 (green on the 



The Effect of Particle Noise in N-body Simulations of Gravitational Lensing 13 



left), 2.9+0.7 (top middle, green+blue) and 8.4 (top right, 
light green) x lO^M©. The influence of the two more massive 
subhalos in the inset with masses of 6.0 and 8.4 xlO^M©, 
exceeds the 3a noise of the critical line. It is evident from 
the critical line of the N-body cluster without any subhalos, 
that the curve still shows a lot of irregularities. Those wig- 
gles are all a consequence of the particle noise due to the 
discrete N-body representation with finite size particles. 

In order to study the infiuence of substructures of dif- 
ferent sizes we could rotate the cluster, however, this would 
also change the overall shape of the critical lines and make 
the comparison between different substructure diflficult to 
quantify. A better approach is then to artificially place ad- 
ditional substructures on top of the critical lines. In the 
bottom panel of Fig. 9 we quantify the effect of artificial 
substructures. We add the substructure particles of 11 sub- 
structures randomly chosen with decreasing lens strength on 
top of the particle distribution of the original N-body sim- 
ulation at ( — 15.6^^—2.6^^). For each substructure we calcu- 
late the new critical lines and compare the deviation from 
the original critical line with the 3a noise on the inverse of 
the magnification in the lower panel of Fig. 9. Substruc- 
tures smaller than /3 ^ 0.1^^ are within the white band in 
the lower panel of Fig. 9. This limit corresponds to resolved 
minimum average subhalo mass of Mss -4.3x lO^M©. Any 
substructure bigger than this will cause a 'wiggle' in the crit- 
ical line that is stronger than those caused by the numerical 
3a particle noise of the simulation. 

6.3 Substructure Deflection Angles 

In this section, we compare the defiection angles of the sub- 
structures ass with the noise on the defiection angles due to 
the particle noise in the N-body simulation, a^- Although 
the defiection angles are not directly observable, we will 
show in this section that they are a good measure of the 
effect of a substructure. 

We have already seen in Sec. 4.5, that the additional 
defiection caused by small-scale fiuctuations due to the par- 
ticle noise can be approximated as a small correction on top 
of the defiection by the numerical cluster. Any additional de- 
fiection Acx will result in a change in the observable surface 
brightness of the image Ad. In a simplified model where 
the source gradient is varying slowly, the greater the ad- 
ditional defiection Aa, the greater the change in the image 
brightness at that point (see Eq. (9)). Therefore we compare 
the additional defiection by a substructure, Aass, with the 
fiuctuations in the defiection angles, a^ caused by the par- 
ticle noise. We could also use an integrated measure over all 
points of the image plane, where the additional defiection by 
the substructure exceeds the magnitude of the fiuctuations 
from the particle noise. Since the results are the same, we 
use here for simplicity only the maximum value of the addi- 
tional substructure defiection, max[Aass {x)]x and compare 
it to acx (x) due to the particle noise at the same point on 
the lens plane. 

Figure 10 shows the maximum of the additional defiec- 
tion for a subsample of the 2597 substructures in the central 
(92^^)^. The subsample is chosen to include the most mas- 
sive substructures with the greatest lens strength. Each sub- 
structure is shown as a red circle for the Level- 2 resolution. 
As a comparison, the subhalos of a second, lower resolution 
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Figure 10. Maximum value of the additional deflection Ao: for 
substructures with positions within (92^^)-^ from the centre for 
the populations Level-2 and Level-4 resolutions. Horizontal lines 
show the Icr (solid lines) and the 3cr (dashed lines) noise at the 
position of the lowest image in the top left panel of Fig. 4 for 
Level-2 (red) and Level-4 (blue). 



of the simulation, Level-4, are also shown as blue squares. 
We use the fully numerical substructures to calculate the 
total mass and half mass radius for the lens strength on 
the X-axis. The y-axis is the maximum of the additional de- 
fiection calculated numerically by solving the Poisson equa- 
tion for each substructure. All subhalos fall almost perfectly 
onto the linear relation y — 0.2272x, therefore the parameter 
Mss/(Scriti?HMR) is sultablc to quantify the strength of the 
numerical substructures and the subsample of 400 subhalos 
is sufficient to quantify the the effects of the subhalo popu- 
lation in Fig. 10. Now, we compare the maximum additional 
defiection caused by the individual substructures with the 
noise on the defiection angles from Sec. 4.4 as follows. We 
artificially place each of the substructures directly behind 
the lowest image at (13^^—4'''). Since the noise on a is a 
very smooth function over the size of a typical image (see 
Fig. 3) we can use the same value a^ for all substructures. 
The noise on the defiection angles can therefore be repre- 
sented by horizontal lines in Figure 10. A more rigorous 
approach would require to use the position of the respective 
maximum of the additional defiection for each substructure 
which varies slightly with the size of the substructure. Using 
the correct appropriate values of ctq-, however, does not alter 
the results. Under the assumption that a greater additional 
defiection by a substructure, Aass, also causes a more sig- 
nificant change in the image brightness, we define as visible 
substructure, those that are above the respective horizontal 
noise levels in Figure 10. 

With this method we are able to constrain the la 'vis- 
ibility' of substructures to a substructure lens strengths of 
0.048'' (1.7 X lO^Mo) for LeveL2 and 0.21'' (1 x lO^^M©) 
for Level-4 (solid lines). The respective 3a limits are 0.14" 
(6.5 X lO^Mo) and 0.62" (3.9 x lO^M©) (dashed lines). This 
simple comparison of the defiections by a substructure with 
the fiuctuations in the defiection from the particle noise pro- 
vides a fast measure of the detect ability of a substructure. 
We will see in the next section that these results based on the 
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Figure 11. Residual R where R = |Ad| — nad > 0, where 
the image brightness difference caused by adding a substructure 
is bigger than the particle noise induced brightness difference in 
the image. Shown is the maximum of R measured over all points 
on the lens plane. Points are for > 2500 substructures from the 
central (92^^)^ of the cluster artificially placed on top of the lowest 
image from the top left panel of Fig. 4. Filled symbols are for a 
fixed source brightness distribution (fs) red circles for la and 
blue squares for 3cr noise, solid lines are a linear fit to the 1 
and 3(7 points. Open symbols are for a reconstructed source (rs) 
where the only assumption is the regularity of the source, black 
squares for la and green circles for 3a. Grey small filled circles are 
the integrated residual over the image plane with a reconstructed 
source (rsi) and la noise. The lower limits of resolved substructure 
lens strengths derived from the linear fits are 0.045''(lcr) and 
0.13''(3cr). 

deflection angles are almost identical to the analysis based 
on differences in the image brightness distribution. 

6.4 Lensed images with Substructure of a 
Gaussian Source 

In this section we compare the effect of substructures within 
the N-body simulation with the particle noise based on the 
surface brightness distribution of the images. The gravita- 
tionally lensed images are the only directly observable quan- 
tities, therefore any observational detection of a substruc- 
ture in the lens will be based on a reconstruction of the ob- 
served image positions or the image brightness distribution. 
If there is no small-scale structure in the lens, it is in theory 
possible to find a source surface brightness distribution that 
fits all of the multiply lensed images. The presence of a sub- 
structure in the lens close to one image of a multiply imaged 
system will result in a local imprint of the substructure on 
the closest image. Therefore, no source brightness distribu- 
tion exists, that, when lensed through a smooth large-scale 
lensing mass distribution, will be able to model all multi- 
ply lensed images simultaneously. Therefore, there has to be 
substructure in the lens. 

Here, we are using a similar method to quantify the ef- 
fect of a substructure. We artificially place a substructure 
on the lens and simulate an image with substructure. We 
then try and fail to reconstruct this lens with substructure 



with a smooth model. In our case, we know the underlying 
smooth model, which is the same lens without the substruc- 
ture. We can therefore calculate the image for the idealised 
smooth model by lensing a source through the lens model 
without the added substructure. The difference between the 
image for a lens with substructure and for a lens without 
substructure is then a measure for the failure of the smooth 
model to reproduce the image with substructure. 

Similarly to Sec. 4.5, we use two different approaches to 
evaluate the influence of substructure on an image bright- 
ness distribution. The first method uses a fixed source sur- 
face brightness distribution to create a reference image do 
using the N-body cluster lens without any substructures 
close to any of the lensed images. This reference image is 
shown in the top left panel of Fig. 4. We then use the 2597 
substructures from the central (92^^)^ as a sample of nu- 
merically simulated substructures and artificially place each 
substructure on top of the cluster directly behind the lowest 
of the three images at (13^^,— 4'^). At this position on the 
lens plane the substructure will have the biggest effect. For 
each substructure we then lens the same source and obtain 
2597 different images di. For each of these images we calcu- 
late the image brightness difference di — do = Adi at each 
point on the lens plane. This is the brightness difference due 
to the artificially added substructure. We then compare this 
brightness difference Adi with the discreteness noise on the 
image ad from Sec. 4.5. We calculate the residual, Ri, as a 
difference in the image brightness due to the substructure 
that is greater than the amplitude of the fluctuations in the 
image surface brightness distribution at that point due to 
the particle noise as 

R^ (x) = \ Ad^ (x) \-nad{x)^0 n = {1, 2, 3} . (24) 

The results for a residual integrated over the lens plane are 
identical to the results from the comparison based on the 
maximum value mdux[Ri {x)]x shown as an example for the 
integrated la residual (rsi) multiplied by a factor of 0.01 to 
enhance the contrast in Fig. 11. The maximum value of the 
residual, Ri, is shown as solid red circles (la) and solid blue 
squares {3a) in Fig. 11 for the Level-2 resolution using a 
fixed source surface brightness distribution (fs). The y-axis 
is the maximum residual Ri and the x-axis is the lensing 
strength of the substructure. The solid lines are linear fits 
that constrain the visible substructures. We consider a sub- 
structure to be visible when Ri > 0. The limits are 0.045^^ 
(Icr), 0.09'' (2cr) and 0.133'' {3a) corresponding to 1.6, 3.8 
and 6 xlO^M© respectively for a fixed source brightness 
distribution. Any substructure with a lens strength greater 
than these lower limits placed on top of the lowest image will 
result in an image brightness difference that is greater than 
the fluctuations due to the particle noise. These limits are 
very close to the values derived from the deflection angles in 
the previous section. 

We have seen in Sec. 4.5 that lensing one fixed source 
surface brightness distribution through different noise reali- 
sations of the lens results in an overestimation of the noise on 
the image by a factor of ^ 2. This definition of the noise on 
the image also includes an artificial relative source-caustic 
shift (see Sec. 4.5 for details), we therefore expect that the 
assumption of a fixed source here also is an over simplifica- 
tion. We therefore use the method described in Sec. 4.5 to 
reconstruct the best possible source (and therefore also the 
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closest image) using the deflection angles with substructure 
to match the same image without substructure. For each of 
the reconstructed substructure images we then calculate the 
residual Eq. (24) from the previous paragraph. The maxi- 
mum of these residuals are shown in Fig. 11 with black open 
squares (la) and green open circles (3cr) for a reconstructed 
source (rs), as well as an integrated residual as an example 
(rsi la). The points no longer follow a linear relation. But 
the lower limits derived with this method of 0.045^^ (la) and 
0.13" (Sa) if we reconstruct the source for each substructure 
image are identical to the ones with a fixed source brightness 
distribution. 

This nontrivial result shows that including the unob- 
servable relative source-caustic shift in the noise increases 
the noise by about a factor of two, see Sec. 4.5. But at the 
same time it also increases the image brightness difference 
due to a substructure in the lens. Therefore the cutoff R = 
which indicates the minimum size of the resolved substruc- 
tures remains unchanged. This is very convenient, since it 
allows us to accurately calculate the limits with much sim- 
pler and faster methods without having to reconstruct the 
source for each of the substructures and bootstrapped parti- 
cle noise realisations of the cluster. The limits for the small- 
est resolved substructures that we found in this section based 
on the image brightness distribution of the lensed images 
are very similar to the limits from simpler lensing proper- 
ties such as the deflection angles in Sec. 6.3 or the surface 
mass density in Sec. 6.1. 

Up to here we have used a Gaussian source surface 
brightness distribution. This simplified source allows for a 
systematic description of the the effect of the particle noise 
on the lensed images, however we expect real source galax- 
ies, especially at z = 2, to be more structured. In the fol- 
lowing, therefore, we simulate sources with different degrees 
of structure. 

6.5 Lensed images with Substructure of a 
Realistic Source 

In order to simulate various degrees of smoothness for the 
source surface brightness distribution, we use the very struc- 
tured source surface brightness distribution of a true galaxy 
as observed with HST. We then smooth this source surface 
brightness distribution with a Gaussian kernel with increas- 
ing sizes denoted as 20, 40 and 100. As an example, the 
source brightness distributions for 100 and 20 are shown in 
the bottom panels of Fig. 7. For each different source we 
recalculate the noise on the image brightness distribution, 
cTcz, as described in Sec. 4.5. We show the residual above 
the particle noise on the image brightness, of Eq. (24) for 
n = 3, in Fig. 12. As the structure in the source is decreased, 
an increasingly smaller number of substructures can be re- 
solved above the noise limit, and the cutoff where R ^ 
shifts to the right. The limit for resolved substructures for 
the smoothest source, 100, is almost identical to the limit 
calculated with a Gaussian source surface brightness distri- 
bution in the previous section. In Fig. 12 we additionally 
show the residual using two different sources with an in- 
creased maximum source surface brightness by factors of 
10 and 100 and the original Gaussian brightness distribu- 
tion. As an extreme and unphysical limit we also show the 
residual for a source surface brightness distribution that is 
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Figure 12. The same as in Fig. 11 but for different source bright- 
ness distributions and 3cr residuals. A very structured source 
brightness distribution is smoothed with a Gaussian kernel of 
sizes 20, 40 and 100, as a comparison, the Gaussian source surface 
brightness distribution from Fig. 11 and a linear source surface 
brightness distribution are also shown. A more structured source 
shifts the cutoff towards smaller substructures, while the overall 
residual is increased. Increasing the brightness of the source by a 
factor of 10 and 100 does not change the smallest resolved sub- 
structures, but increases the residual. For the sake of the clarity 
of the figure, only every 5th substructure is plotted for each line. 



decreasing linearly with increasing distance from the source 
centre. Size and position of the source are kept constant for 
each of the seven curves. The behaviour of the curves can 
be qualitatively understood as follows. 

Changing the brightness of the source while keeping the 
source size constant, increases the gradient of the source. 
Therefore, the same change in the deflection angles, Aa, 
due to either noise or substructure, will result in a greater 
change in image brightness. Ad ^ Aa • Vs, see also Eq. (9). 
This effect, however, affects the calculation of the noise on 
the images, 3crd, and the image brightness difference for 
each substructure, Adfc, in the same way. Therefore even 
though the residual Eq. (24) is increased, the lower limit 
of resolved substructures is unchanged. To understand the 
behaviour of the curves with increasing source structure, 
we consider a small substructure that is barely not resolved 
with a Gaussian source brightness distribution, for example 
a substructure with a lens strength of 0.12^^ which corre- 
sponds to a substructure mass of - 5.4 X lO^M©. bmce 
the residual of this substructure is i?max ^ 0, the image 
brightness difference at any point due to this substructure, 
Adi, is comparable to the noise on the image brightness, 
3crd, which we calculated from the different particle noise 
realisations of the cluster, Gd = 1/(A^ — 1)^^ /\dk. If we 
now increase and change the form of the source gradient 
with a more structured source on scales of Aa^, the same 
change in the deflection angles due to the substructure, Aa^, 
will result in an increased change in the image brightness, 
Adi. In contrast to the increased source brightness, which 
only changes the amplitude of the source gradient, a more 
structured source additionally changes the shape of the gra- 
dient on scales smaller than Aa^. Therefore and because 
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Figure 13. Required number of particles in a N-body simula- 
tion targeted to resolve substructures above 3cr of the noise with 
a given lens strength with strongly gravitational lensed images. 
The second x-axis converts the lens strength of the substructure to 
the average mass of a simulated substructure. The horizontal lines 
show three resolutions of the Phoenix simulations and their re- 
spective lower 3fj hmits of ^3™^^ - 0.046'' or M^'"^ - 1.7x lO^M© 
for Level-1 and ~ 0.135'' or M™^^ 5.8 x lO^M© for Level-2 
and 13^''^ - 0.6" or M^'"^ - 3.8 x IO^OM© for Level-4. 



the noise ad is a nonlinear combination of the effects of a 
field of positive and negative noise substructures it behaves 
differently than a more coherent single substructure. As we 
can see from Fig. 12 the net effect is more prominent for 
a single substructure. An increased source structure there- 
fore allows the detection of smaller substructures above the 
particle noise limit. 

Fig. 12 shows that a Gaussian source surface brightness 
distribution is a bad choice in terms of the lower limit of re- 
solved substructures. The results are almost identical to the 
worse case, a linear source brightness distribution The truly 
worst choice, however, would be a constant source brightness 
which is completely indifferent to changes in the deflection 
angles. A linear source gradient includes both, an increased 
gradient and at the same time less structure in the source 
brightness distribution with respect to a Gaussian source 
surface brightness distribution. Therefore, the overall resid- 
ual is slightly higher and the cutoff is shifted to greater sub- 
structure lens strengths. The linear source surface brightness 
distribution is the extreme limit, a constant source gradient 
in Eq. (9). Therefore, the cutoff at /3 = 0.14'' (6.5 x lO^M©) 
is identical to the lower limit derived from the substructure 
deflection angles in Sec. 6.3. 



7 SCALING OF THE RESOLUTION LIMIT 

In the previous sections we derived different methods to 
compare the effect of a substructure on different lensing 
properties with the effect caused by the discrete representa- 
tion of the N-body cluster with particles and we found lower 
limits on the substructure lens strength. Substructures with 
a weaker lens strength do not affect the lensing properties 



strong enough in order to be detectable above the particle 
noise limit. 

This allows us to answer a very interesting question. If 
we plan to simulate gravitational lensing with a numerically 
simulated N-body mass distribution, what resolution do we 
need, in order to 'resolve' substructures of a given size. 

The most advanced lower limits derived from the image 
brightness distribution with a reconstructing source method 
are very close to the simplified limits predicted from deflec- 
tion angle differences, or even the simple comparison based 
on the surface mass density fluctuations in Sec. 6.1. We al- 
ready derived scaling relations for the noise on the surface 
mass density in Sec. 5, therefore, we can now combine those 
results. If we approximate the noise as positive and negative 
Gaussian fluctuations, we can write the lens strength of one 
noise fluctuation as /3 = 7r/cr,^/(2 In 2), see Sec. 6.1. Here, / 
is the size of the smoothing kernel that is used to smooth 
the N-body particles on the lens plane, and cr^t is the noise 
on the surface mass density. We can use this approximation 
to estimate the scaling of the equivalent lens strength of the 
noise fluctuations with A^part and A^ngb- Using Eqs. (19) and 
(20), the relation scales as 
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(25) 



The limit for resolved substructures thus is independent of 
A^ngb- This approximation holds, as long as the increase in 
smoothing does not smooth out any small-scale substruc- 
tures and as long as the surface mass density remains rea- 
sonably smooth. To test Eq. (25), we repeated the calcula- 
tions in Sec. 6.4 which were done with a smoothing kernel 
with A^ngb = 64, this time with a reduced number of neigh- 
bours A^ngb = 8- This yields the same lower limit for the 
resolved substructures, however the noise on the individual 
lensing properties is substantially more prominent for a less 
smoothed N-body particle distribution. In Sec. 4 we quanti- 
fied the noise individually for each lensing quantity. For ex- 
ample the critical lines and therefore the caustic will exhibit 
a great number of higher order singularities and swallowtails 
if the smoothing of the N-body particles is reduced. To sim- 
ulate gravitational lensing, we have to smooth the N-body 
particles in order to obtain a reasonably smooth simulated 
image that is free from too much artificial numerical sub- 
structure. 

From Sec. 6.4 we know that the Level-2 resolution re- 
solves substructures as small as /3 ~ 0.135" . Inserting this 
result in Eq. (25), we can quantify the resolution require- 
ment for a N-body simulation of gravitational lensing as 
a function of the smallest resolved substructures, which is 
shown in Fig. 13. For a given substructure that we want to 
resolve, we can estimate the size of the simulation we need 
to reduce the noise enough in order to see an effect of the 
chosen substructure above 3a of the particle noise. Figure 13 
is calculated from the noise in the surface mass density in 
Sec. 6.1, but the limits for the deflection angles in Sec. 6.3 
and the lensed images with a linear source surface brightness 



The Effect of Particle Noise in N-body Simulations of Gravitational Lensing 17 



distribution in Sec. 6.4 are identical. We have seen in Fig. 12 
that the true lower resolution limits depend on the detailed 
source structure. However, we also know from Fig. 12 that 
a linear source surface brightness distribution is the worst 
case scenario. Therefore, these limits are valid, independent 
of the source used to simulate the gravitational lensing. 



8 SUMMARY 

In the first part of this paper we investigated the effect of the 
discrete representation of a N-body simulation with particles 
on the simulation of gravitational lensing. We have used the 
currently highest-resolution simulations, the Phoenix simu- 
lations for our numerical lensing simulation. With the reso- 
lution of the Level-2 simulation we found the noise on the 
projected surface mass density to be 1 — 2% and smaller than 
5% even in high density regions at the centre of the cluster. 
We then used the noise on the inverse of the magnification 
to quantify the irregularities of the critical line. The noise 
in units of the inverse of the magnification increases with 
magnification and reaches its maximum at the critical lines 
where the magnification diverges. Due to the nonlinearity 
of the lensing equation the noise plays a significant role in 
these high-magnification regions of strongly lensed images. 
In the surface brightness distribution of the lensed images, 
the particle noise causes fiuctuations of the order of 10%. 
However, the assumption of a static source surface bright- 
ness distribution for the calculation of the noise fiuctuations 
includes some unobservable effects such as a shift of the de- 
fiection angles relative to the source. Therefore, we also used 
a Bayesian source reconstruction argument in order to prop- 
erly quantify the noise on the multiply lensed images. With 
this second method, the particle noise still has a consider- 
able effect on the morphology of strongly magnified images. 
For the Level-2 resolution of the Phoenix simulations and a 
typical three image highly magnified giant arc, the noise on 
the image brightness is ^ 5%. In Section 5 we derived useful 
scaling relations for the particle noise on the surface mass 
density with the number of particles in the simulation and 
the number of smoothing neighbours. 

In the second part of the paper we compared the in- 
fiuence of physical substructures in the N-body simulation 
with the particle noise we derived in the first part of the pa- 
per. We compared the projected surface mass density and 
the deviations of the critical lines. From the comparison of 
the lensed images we found that for substructures with a 
lens strength smaller than 0.13^^ there is no measurable ef- 
fect in a simulation comparable to our Level-2 resolution 
above 3cr of the effect of the numerical particle noise. A typ- 
ical substructure with a lens strength of 0.13^^ has a mass of 
~ 6 X 10^ M© and therefore consists of ^ 10^ particles with 
a mass of 6 x 10^, see Table 1. These results were calculated 
with a fully adaptive source to avoid unobservable effects 
such as a shift in the effective source position. This measure 
of the importance of a substructure is motivated by obser- 
vational reconstruction techniques. We found, that the re- 
sults with the fully Bayesian source reconstruction measure 
are comparable to the much simpler results obtained from 
a non-adaptive, fixed source brightness distribution. In fact, 
the simpler comparisons based on the additional defiection 



caused by a substructure or the effective lens strength of the 
noise mass density fiuctuations yield comparable results. 

Therefore, finally, we combined in Sec. 7 the analysis of 
the scaling of the noise with the particle number of the sim- 
ulation and the number of smoothing neighbours from the 
first part of the paper with the investigations of the effect 
of the simulated substructures from the second part of the 
paper. This allowed us to quantify the required resolution 
of a numerical N-body simulation if we want to detect sub- 
structures of a given size in our simulation of gravitational 
lensing. 



APPENDIX A: VARIANCE OF A SMOOTHED 
PARTICLE DISTRIBUTION 

We calculate the variance of a smoothed distribution a\ (cc) 
from a known variance at each point Uy {x) under the as- 
sumption the latter is a slowly varying function with x. The 
smoothing acts as a convolution with a normalized kernel 
function, which depends on one parameter, the smoothing 
length / (cc) = Ix. For a Gaussian smoothing kernel the value 
of the property we are interested in, X, at the point x can 
therefore be written as 



X{x, h) = j d^x' Y [x') I ^ exp 



{x-x'f 

(Al) 

For random and uncorrelated variables Y (x^) we therefore 
get the variance 



CTxixJx) 



j2 / 2 / / 

d X ayix 



1 



exp 



2/2 



If we now assume that ayix' ^Ix') is 



(A2) 

. slowly varying func- 



tion with position x^ or more precise we assume ay to be 
constant with respect to changing x' over sizes where the 
squared Gaussian kernel is non- vanishing, which is approxi- 
mately true for |cc — cc'l > 3/, we can simplify Eq. (A2) 



2 / 

cfy(x 



1 



(Ty{x, 



47r/2 



^CFy{x,Ix) 



exp 



{x-x'f 



2/2 



(A3) 



The exact value of the constant c, will depend on the form 
of the smoothing kernel, here for a Gaussian kernel c— 1/4. 
As described in Sec. 4.2 there are other smoothing kernels 
with a similar shape, for example the polynomial kernel 
in Eq. (4) with c = 1030/343. We are using the kernel 
Wi{r) = 3/(7r/f) (1 - r'^/l'iY for r <h with c = 9/5. When 
comparing the different numerical values, we have to keep in 
mind that the characteristic length of the Gaussian kernel 
should be decreased by y^TOS/lTSO ~ 0.303 in order for a 
particle to cover approximately the same area compared to 
the other two kernels. 
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